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ABSTRACT 

Strategies  used  to  verify  and  validate  computational  fluid  dynamics  (CFD)  calculations  are  described  via  case 
studies  of  realistic  flow  simulations,  each  representing  a  complex  flow  physics  and  complex  geometry.  Critical 
areas  of  importance  to  validation  of  a  calculation  are  pointed  out  through  various  high  fidelity  physics-  and 
engineering-based  simulations.  These  areas  include  the  physical  model,  conceptual  model,  boundary 
conditions,  initial  conditions,  geometry,  grid  density  and  distribution,  turbulence  model,  and  numerical 
dissipation.  Appropriate  selection  and  exercise  of  the  above  items  depend  upon  thorough  understanding  of  the 
physics  of  the  problem  and  require  considerable  experience  on  addressing  them  utilizing  CFD  codes.  The 
cases  presented  include  the  most  frequently  encountered  flow  features  of  separation,  swirl  and  rotation  as 
observed  in  engineering  applications  such  as  aircraft  gas  combustors  and  turbomachinery.  Each  simulation 
considers  the  salient  physical  features  involved  and  resolves  them  to  the  level  required  by  the  purposes  for 
which  they  are  being  used.  In  many  of  the  simulations  presented,  unstructured  grids  and  parallel  computing 
are  used  to  minimize  the  overall  time  needed  to  achieve  a  numerical  solution.  In  these  cases,  the  numerical 
scheme  incorporated  allows  use  of  a  large  number  (thousands)  of  processors  in  parallel,  to  shorten  the 
solution  time  and  to  provide  speed-ups  that  do  not  deteriorate  with  the  addition  of  more  processors. 
Verification  studies  compare  the  exact  analytical  solutions  with  those  obtained  using  various  numerical 
schemes.  In  these  studies  preservation  of  a  vortex  convecting  through  a  calculation  domain  is  used  to  assess 
the  numerical  accuracy  of  time-dependent  numerical  schemes.  Magnitude  and  distribution  changes  of  field 
variables,  primarily  pressure,  are  quantified  to  obtain  accurate  assessment  of  numerical  errors  with  first- 
order  time-accurate  central-difference,  second-order  time-accurate  iterative  implicit,  and  fifth-order  accurate 
upwind-biased  schemes.  A  scheme  which  minimizes  numerical  error  involves  preservation  of  both  the  vortex 
strength  and  the  vortex  structure. 


1.0  INTRODUCTION 

A  considerable  amount  of  work  has  gone  into  developing  metrics  for  validation  and  verification  of  the 
computational  codes  that  are  developed  and  used  for  modelling  and  simulation  of  fluid  flow  [1-12].  American 
Institute  of  Aeronautics  and  Astronautics  (AIAA),  and  American  Society  of  Mechanical  Engineers  (ASME) 
have  declared  policy  statements  and  guidelines  for  the  verification  and  validation  of  computational  fluid 
dynamics  simulations  [13-16].  These  metrics  and  guidelines  cover  issues  such  as  assessment  for  the  iterative 
convergence,  spatial  grid  convergence,  temporal  convergence,  comparison  of  the  CFD  results  to  experimental 
data,  uncertainty  analysis,  and  assessment  of  the  computer  programs  or  codes  to  check  for  coding  errors,  to 
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mention  a  few.  Fulfilling  all  the  above  mentioned  metrics  and  following  all  the  guidelines  still,  does  not 
provide  a  guarantee  that  a  code  is  truly  verified  and  validated.  As  Roche  [13]  mentioned  it  may  never  be 
appropriate  to  say  that  a  code  itself  is  validated.  The  important  factor  is  the  experience  and  the  ability  of  the 
user  in  understanding  the  underlying  physics  that  govern  a  given  flowfield  and  his  subsequent  evaluation  on 
whether  a  code  can  reliably  model  that  physics.  Otherwise,  the  code  needs  to  be  validated  for  the  given 
problem  in  hand. 

This  paper  describes  a  number  of  numerical  simulations  for  realistic  cases  where  an  assessment  of  the 
underlying  physics  of  the  flow  is  accompanied  with  a  methodical  approach  on  examining  whether  the  code 
can  simulate  those  physics  or  not,  prior  to  applying  the  code  for  the  eventual  simulation. 

2.0  VORTEX  PRESERVATION  TEST 

Vortex  preservation  test  can  be  used  to  validate  the  capability  of  a  code  for  time  accurate  preservation  of  the 
vortical  structures  in  the  flow,  and  the  strength  of  those  vortical  structures  as  they  are  convected  by  the  flow 
downstream.  Vortical  structures  prevail  in  a  number  of  realistic  flows.  The  performance  of  various  devices 
and  vehicles  are  dramatically  affected  by  these  vortical  structures.  It  is  therefore  crucial  that  the  codes  used  in 
these  applications,  where  vortical  structures  play  a  role,  are  capable  of  maintaining  these  structures.  Examples 
of  the  flows  where  preserving  vortical  structures  is  essential  for  meaningful  results  are  flowfields  of 
hurricanes,  tornados,  blade  vortex  interaction  (BVI),  manoeuvring  submarines,  and  reacting  flows  in 
combustion  chambers. 

Preservation  of  a  vortex  convecting  through  a  computational  domain  is  used  here  to  assess  the  numerical 
accuracy  of  time-dependant  Navier-Stokes  numerical  schemes.  Magnitude  and  distribution  changes  of  field 
variables,  primarily  pressure,  are  quantified  to  obtain  accurate  assessment  of  numerical  errors  with  first-order 
time-accurate  central-difference,  second-order  time-accurate  iterative  implicit,  and  fifth-order  accurate 
upwind-biased  schemes.  During  convection  of  a  Lamb-type  vortex  [17]  in  a  freestream,  pressure  is  at 
minimum  in  the  vortex  core  and  increases  concentrically  and  asymptotically  to  the  freestream  value  with 
radial  distance  from  the  vortex  center.  A  numerically  dissipative  scheme  cannot  maintain  this  pressure 
minimum  at  the  vortex  center,  and  core  pressure  increases  as  the  vortex  convects  downstream.  In  inviscid 
flow  calculations  such  as  those  solving  Euler  equations,  the  pressure  distribution  should  remain  unaltered  by 
convection  but  numerically  dissipative  schemes  result  in  a  predicted  pressure  which  increases  with 
downstream  distance  caused  by  the  numerically  induced  vortex  decay.  The  concentric  structure  of  the  vortex 
should  also  remain  unaltered  with  downstream  convection,  however  certain  numerical  schemes  cannot 
maintain  concentricity  of  isobars  around  the  vortex,  instead  patterns  are  distorted  and  characterized  by 
solution  oscillation  waviness.  A  scheme  which  minimizes  numerical  error  involves  preservation  of  both  level 
and  distribution  of,  for  example  the  pressure  minimum  and  concentricity  of  the  vortex,  respectively. 
Preservation  of  concentricity  alone  may  be  insufficient  and  misleading. 

A  series  of  calculations  using  a  conventional  first-order-accurate,  central-difference  scheme  and  a  second 
order  unsteady  iterative  implicit  scheme  are  performed  by  solving  the  time  dependent  Euler/Navier-Stokes 
equations  [18].  Variation  of  vortex  core  pressure  against  number  of  core  radii  travelled  by  the  vortex  is 
monitored  to  measure  the  time  and  spatial  accuracy  of  the  Navier-Stokes  code.  These  calculations  are 
performed  to  assess  the  vortex  preservation  capabilities  of  the  numerical  procedure  and  to  indicate  areas 
where  refinement  would  be  required.  Based  on  the  results  of  these  calculations,  the  code  is  further  developed 
to  provide  the  ability  to  perform  practical  studies  of  airfoil- vortex  interaction  [19].  Details  of  the  numerical 
scheme,  governing  equations  and  the  vortex  definition  are  given  in  [18]. 
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2.1  COMPUTATIONAL  GRID 

The  computational  domain  used  for  the  convection  of  the  vortex  is  shown  in  Figure  1.  The  grid  is  non- 
uniform,  with  equal  spacing  (Ax  =  Ay  =  1/8)  used  in  the  central  region  containing  the  vortex  path  and  a 
stretched  grid  is  used  to  extend  the  physical  domain  to  the  far  field.  The  boundaries  of  the  grid  are  245  x  200 
radii  in  length  to  width,  the  core  radius  of  the  vortex  being  1.0.  The  boundaries  of  the  equally  spaced  grid 
region,  the  overall  boundaries  and  the  vortex  path  are  shown  in  Figure  1.  The  boundaries  of  the  grid  are 
considerably  far  from  the  vortex  path  to  eliminate  any  boundary  effects  on  the  vortex. 
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Figure  1:  Physical  domain  and  the  vortex  path 


2.2  FIRST  ORDER  IN  TIME,  CENTRAL  DIFFERENCE  SCHEME 

A  calculation  is  performed  to  determine  the  flow  characteristics  associated  with  a  Lamb-type  vortex 
convection  in  a  freestream,  using  the  central  difference  scheme  with  first  order  accuracy  in  time,  without  any 
inner  iteration  at  each  time  step.  In  these  calculations,  the  reference  length  is  the  vortex  core  radius,  the 
reference  flow  conditions  are  the  free  stream  conditions  with  the  Mach  number  M  =  0.536.  Time  t,  is  made 
dimensionless  by  freestream  velocity  and  the  vortex  core  radius,  e.g.  an  increment.  At  =  1.0,  is  the  time 
required  for  a  particle  at  freestream  velocity  to  travel  one  vortex  core  radius.  Curve  A  in  Figure  2  shows  the 
variation  of  the  vortex  core  pressure  with  the  distance  travelled  for  this  first-order  scheme,  second-order 
iterative  implicit,  and  fifth-order  upwind-biased  scheme.  As  shown  by  the  curve  A  in  Figure  2,  the  predicted 
core  pressure  is  increased  drastically  compared  to  its  initial  value.  It  is  evident  that  the  first-order  time 
accurate  scheme  is  very  dissipative;  an  improvement  is  needed  for  vortex  preservation.  This  is  achieved  by 
using  the  second-order  unsteady  iterative  implicit  scheme  described  in  the  next  section.  However,  the  results 
of  this  calculation  are  consistent  with  those  obtained  by  Rai  [20]  who  used  the  first-order  time  accurate  Beam 
and  Warming  scheme.  It  should  be  noted  that  although  the  vortex  has  lost  a  significant  amount  of  its  strength, 
the  vortex  shape  is  very  well  preserved  after  45  radii  of  travel. 
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Fiyiire  3:  Pressure  ceiitoiirs  at  Initializatien  (top) 
Figure  2:  Vortex  core  pressure  decay  rate  after  vortex  has  travelled  45  radii  (hottom) 

three  diftereiit  schemes.  2iid-order  iterative  implicit  scheme 


2.3  SECOND-ORDER  TIME  ACCURATE  ITERATIVE  IMPLICIT  SCHEME 

The  basic  scheme  used  is  a  Linearized  Block  implicit  ADI  procedure  of  Briley  and  McDonald  [21].  The 
splitting  error  and  the  linearization  error  associated  with  this  basic  scheme  are  removed  by  introducing  an 
inner  ADI  iterative  procedure  at  each  time  step.  The  temporal  accuracy  is  increased  to  second-order  by  using 
three-point  backward  time  differencing.  On  convergence  of  the  inner  iteration,  the  scheme  becomes  a  fully 
implicit  nonlinear  backward  time  difference  scheme.  A  more  detailed  discussion  of  these  improvements  is 
given  by  Rai  [20]  and  such  an  iterative  implicit  three-time  level  ADI  scheme  is  used  in  the  present  study.  In 
addition,  three  point  central  differences  are  used  to  approximate  the  spatial  derivatives.  The  spatial  accuracy  is 
second-order  except  for  the  use  of  numerical  dissipation  discussed  subsequently. 

When  calculating  high  Reynolds  number  using  either  the  Euler  equations  or  the  Navier-Stokes  equations, 
using  centred  spatial  differencing,  some  artificial  dissipation  is  usually  needed  to  maintain  numerical  stability 
and  to  suppress  spurious  oscillations  in  the  numerical  results.  The  approach  used  in  the  present  effort  is  based 
upon  the  use  of  a  second-order  anisotropic  artificial  dissipation  term.  Introduction  of  the  second-order  terms 
for  artificial  dissipation  formally  reduces  the  scheme  to  first-order.  However,  the  added  second-order  artificial 
dissipation  term  is  preceded  by  an  adjustable  parameter  which  can  be  reduced  so  as  to  progressively  reduce 
the  effect  of  this  term.  The  parameter,  termed  AVISC  in  the  figures  which  will  be  presented,  is  essentially 
equivalent  to  a  coordinate  dependent  inverse  cell  Reynolds  number  so  that  a  specification  of  AVISC  of  0.05 
limits  the  maximum  cell  Reynolds  number  to  20.  Obviously,  artificial  dissipation  is  a  source  of  false  diffusion 
and  distortion  of  a  vortex  convected  in  a  finite-difference  grid.  In  the  present  work,  the  effects  of  artificial 
dissipation  on  the  vortex  structure  is  minimized  by  specifying  a  small  value  of  the  adjustable  parameter  which 
controls  the  amount  of  added  dissipation,  and  this  value  is  determined  from  a  separate  set  of  calculations  in 
which  the  effects  of  this  parameter’s  magnitude  on  the  preservation  of  the  free  vortex  convected  over  a  long 
distance  are  examined.  The  values  used  for  this  parameter  have  some  effect  on  the  results,  as  will  be 
demonstrated.  In  some  cases  a  value  of  zero  could  be  used  but  it  was  not  determined  that  this  was  true  in  all 
cases.  Consequently  for  safety  a  small  valve  of  AVISC  was  used  routinely.  The  variation  of  vortex  core 
pressure  with  number  of  core  radii  is  shown  by  curve  C  in  Figure  2,  which  indicates  a  very  stable  and  accurate 
solution.  There  was  no  significant  distortion  of  the  vortex  during  this  travel.  Contours  of  pressure  at 
initialization  and  after  45  radii  of  the  vortex  travel  are  also  shown  in  Figure  3. 
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2.4  EFFECTS  OF  SPATIAL  SPACING  ON  THE  SOLUTION  ACCURACY 


Two  calculations  were  performed  on  grids  with  different  spacing  to  study  the  effects  of  grid  spacing  on  the 
vortex  preservation.  The  first  grid  has  the  spacing  of  Ax  =  Ay  =  %  ,  and  the  second  grid  is  a  finer  grid  with  the 
spacing  of  Ax  =  Ay  =  1/8.  Variation  of  vortex  core  pressure  versus  number  of  core  radii  traveled  for  both  grids 
is  shown  in  Figure  4.  It  can  be  seen  that  numerical  solution  for  the  coarser  grid  is  oscillating,  particularly 
toward  the  latter  parts  of  the  vortex  travel.  Contour  plots  of  pressure  after  45  radii  (not  shown)  showed  a  badly 
deformed  vortex.  Variation  of  vortex  core  pressure  with  number  of  core  radii  for  the  finer  grid  shown  by  the 
curve  B  in  Figure  4,  indicates  a  very  stable  and  accurate  solution.  In  contrast  to  the  coarser  mesh  case,  this 
case  showed  no  significant  distortion  of  the  vortex  during  its  travel. 
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Figure  4:  Effects  of  the  grid  spacing 
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Figure  5:  Effects  of  the  Numericai  Dissipation 


2.5  EFFECTS  OF  NUMERICAL  DIFFSUSION 

Use  of  artificial  diffusion  enhances  the  stability  and  convergence  properties  of  the  numerical  solution 
procedures.  Such  artificial  diffusion  could  be  added  via  the  spatial  differencing  formulation  (e.g.  one-sided 
difference  approximations  for  first  derivatives)  or  by  explicitly  adding  an  additional  diffusion  term.  In  the 
numerical  scheme  used  in  this  study  the  latter  approach  was  adopted,  since  when  an  additional  term  is 
explicitly  added,  the  physical  approximation  being  made  is  usually  clearer  than  when  dissipative  mechanisms 
are  contained  within  numerical  truncation  errors,  and  further,  explicit  addition  of  an  artificial  viscosity  term 
allows  greater  control  over  the  amount  of  non-physical  numerical  dissipation  being  added.  Obviously,  the 
most  desirable  technique  would  add  only  enough  numerical  dissipation  to  suppress  oscillations  without 
deteriorating  solution  accuracy.  Four  cases,  including  zero  artificial  diffusion,  were  run  to  assess  the  effect  of 
the  numerical  diffusion  on  the  unsteady  vortex  flow  and  to  help  choose  values  which  would  preserve  accuracy 
and  yet  suppress  the  oscillation.  Variation  of  the  vortex  core  pressure  versus  number  of  core  radii  travelled  by 
the  vortex  for  different  values  of  the  artificial  viscosity  is  given  in  Figure  5.  The  parameter  “AVISC” 
indicated  in  the  figure  is  a  measure  of  the  amount  of  artificial  viscosity  added.  The  larger  the  AVISC  value, 
the  larger  the  added  artificial  viscosity.  A  value  of  0.5  corresponds  to  a  first  order  upwind  difference  level  of 
numerical  dissipation.  The  values  used  are  0.0,  0.001,  0.005,  and  0.050.  It  is  clear  from  these  curves  that  the 
higher  the  value  of  the  artificial  viscosity  the  more  dissipative  the  calculation.  It  should  be  noted  that  in  the 
present  calculation  the  solution  remains  stable  without  any  added  artificial  viscosity  throughout  this 
calculation.  Contours  of  pressure  and  vorticity  magnitudes  after  45  radii  of  travel  for  calculations  with 
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different  values  of  artificial  viscosity  indicated  that  the  vortex  shape  is  very  well  preserved  after  45  radii  of 
travel,  despite  the  fact  that  for  higher  values  of  artificial  viscosity  the  vortex  has  lost  more  of  its  strength. 


2.6  DIAGONAL  CONVECTION 

This  case  is  conducted  to  evaluate  the  capability  of  the  code  in  preserving  the  vortical  structures  travelling 
diagonal  to  the  grid  lines  and  in  a  direction  not  aligned  with  a  coordinate  direction.  With  the  capability  to 
convect  a  vortex  in  any  direction,  one  could  study  the  interaction  of  a  vortex  flowing  into  a  blade  in  any 
arbitrary  direction. 

The  computation  was  done  using  the  second  order  accurate  iterative  implicit  scheme.  A  lamb  type  vortex  as 
described  2.0  is  initialized  in  a  bottom  corner  of  a  rectangular  domain  and  is  convected  by  the  flow  diagonally 
to  the  opposing  corner.  Contours  of  the  static  pressure  at  the  initial  condition  are  shown  in  Figure  6,  with  the 
inflow  angle  at  a  45  degree  angle  to  the  x-coordinate  direction.  Contours  of  the  pressure  after  the  vortex  has 
travelled  20  radii  and  45  radii  are  also  shown  in  this  Figure.  Grid  size  for  this  calculation  was  337  x  337.  It  is 
noticed  that  the  vortex  has  maintained  its  structure.  Furthermore,  monitoring  the  vortex  core  pressure  along 
its  path  was  similar  to  the  curve  C,  shown  in  Figure  2,  indicating  that  the  vortex  has  preserved  its  strength  in  a 
diagonal  to  gird  lines  path,  under  this  numerical  scheme. 


3.0  BLADE  VORTEX  INTERACTIONS  (BVI) 

The  interaction  of  concentrated  vortices  with  blades  induces  unsteady  aerodynamic  loading  responsible  for 
blade  vibrations,  aeroelastic  instabilities,  and  impulsive  noise.  The  effects  of  blade-vortex  interaction  (BVI) 
are  especially  significant  in  the  transonic  flow  regime,  in  which  the  strength  and  position  of  the  shock  waves 
are  sensitive  to  small  changes  in  the  flow  parameters.  At  the  present  time,  a  key  problem  in  computing  flows 
containing  concentrated  vortices  is  the  ability  to  preserve  and  convect  these  vortices  in  a  finite  difference  or 
finite-volume  grid  without  false  numerical  diffusion  due  to  truncation  error,  artificial  dissipation  and 
turbulence  modelling.  In  this  work,  the  ensemble-averaged,  time-dependent  Navier-Stokes  equations  are 
solved  on  a  body-fitted  grid  around  a  NACAO012  airfoil  to  study  strong  interaction  of  a  vortex  with  a 
stationary  blade.  The  Navier-Stokes  equations  are  solved  by  using  an  iterative  implicit  finite  difference 
scheme  with  second  order  spatial  and  temporal  accuracies.  Furthermore,  simple  vortex  preservation 
techniques  discussed  earlier  are  used  to  minimize  the  amount  of  spurious  numerical  dissipation  and  eddy 
viscosity  caused  by  the  presence  of  the  vortex  during  its  convective  motion  towards  the  leading  edge  of  the 
blade.  If  a  numerical  procedure  does  not  preserve  an  isolated  vortex,  this  same  procedure  can  not  be  used  for 
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the  blade  vortex  interaction  study  since  non-physical  dissipation  will  occur  prior  to  the  interaction  process  and 
the  process  will  be  based  upon  a  vortex  of  artificially  low  strength  and  large  extent. 


3.1  FLOW  PARAMETERS 

The  reference  length  is  the  chord  length  of  the  blade  and  the  reference  flow  conditions  are  the  free  stream 
condition  with  Moo  =  0.8  and  Re  =  1.0  x  10^.  The  background  flow  is  a  steady  transonic  flow  with  shock 
waves  standing  in  the  middle  of  the  blade.  Furthermore,  the  flow  is  symmetric  about  the  chord  line;  hence,  the 
lift  coefficient  (CL)  is  zero.  The  surface  pressure  distribution  of  this  background  flow  is  shown  in  Figure  7. 
The  dimensionless  strength  and  core  radius  of  the  vortex  are  -1.6  and  0.2,  respectively,  where  the  minus  sign 
indicates  that  the  vortex  has  a  clockwise  sense.  The  initial  location  of  the  vortex  center  is  at  a  point  5  chords 

upstream  of  the  airfoil  leading  edge  (x^  =  -5.0)  and  0.26  chords  below  (Y^  =  -0.26).  The  calculation  is  carried 
out  from  t  =  0  to  t  =  8  with  constant  time  step  At  =  0.005.  It  is  noted  here  that  the  vortex  core  arrives  at  the 
blade  leading  edge  when  t  =  4.95,  which  indicates  an  average  core  velocity  of  0.99  Voo. 

The  total  number  of  grid  points  used  is  144  x  1 18.  The  inflow  boundary  is  located  at  7  chords  from  the  blade 
leading  edge  while  the  outflow  boundary  is  located  at  5  chords  from  the  blade  trailing  edge.  Based  upon  the 
isolated  vortex  study  the  distance  between  the  top  boundary  and  the  chord  line  of  the  NACAO012  airfoil  was 
set  at  5  chord  lengths.  The  geometric  configuration  is  symmetric  about  the  chord  line.  Along  the  inflow 
boundary,  the  total  pressure,  the  total  temperature  and  the  inflow  angle  are  specified.  Along  the  outflow,  top 
and  bottom  boundaries,  the  static  pressure  is  specified;  the  velocity  and  the  total  temperature  are  obtained  by 
extrapolation.  On  the  blade  surface,  non-slip  conditions  are  imposed.  The  density  is  obtained  by  solving  the 
continuity  equation  and  the  surface  temperature  is  specified  as  the  constant,  free  stream  total  temperature. 

3.2  BVI  SIMULATION  RESULTS 

The  interactions  between  the  vortex  and  the  blade  with  a  shock  are  further  elucidated  in  terms  of  the 
instantaneous  static  pressure  distribution  at  several  selected  time  stations.  Figure  6  gives  the  pressure  contours 
over  the  entire  computational  domain  at  t  =  0.  Distribution  of  static  pressure  coefficient  on  the  blade  surface  at 
t  =  0  is  also  shown  in  this  Figure.  This  is  the  starting  flow  field.  As  the  vortex  convects  towards  the  blade,  the 
upper  surface  shock  moves  in  the  upstream  direction,  its  strength  is  decreasing  and  the  extent  of  the  associated 
supersonic  pocket  also  is  reducing.  On  the  other  hand,  the  lower  surface  shock  moves  in  the  downstream 
direction  with  increased  strength.  In  addition  to  the  motion  of  the  shock  waves,  pressure  difference  between 
the  upper  and  lower  surfaces  start  to  build  up.  These  generic  features  are  illustrated  in  Figures  8  and  9.  The 
upper  surface  supersonic  pocket  practically  has  disappeared.  The  lower  surface  shock  wave  becomes  stronger 
and  is  located  in  a  further  downstream  position;  at  the  shock's  root  the  flow  shows  signs  of  separation.  The 
emission  of  a  high  pressure  pulse  from  the  upper  surface  of  the  leading  edge  is  evident  from  Figure  9;  this 
high  pressure  pulse  then  propagates  upstream  in  a  domain  including  the  frontal  region  of  the  entire  leading 
edge.  Between  this  frontal  high  pressure  region  and  the  lower  surface  shock  wave,  a  low  pressure  pulse  is 
propagating  towards  the  lower  outer  boundary.  The  general  features  of  the  flow  at  t  =  6.0  are:  the  existence  of 
a  supersonic  pocket  on  the  lower  surface,  significant  flow  separation  originating  at  the  root  of  the  shock,  the 
appearance  of  vortex  remnants  near  the  blade  trailing  edge,  and  the  development  of  supersonic  flow  on  the 
upper  surface.  The  flow  on  the  lower  surface  does  not  exhibit  any  appreciable  separation  and  is  entirely 
transonic.  Furthermore,  about  70%  of  the  upper  surface  is  covered  by  a  supersonic  pocket,  with  compression 
waves  appearing  near  the  trailing  edge  of  the  blade.  It  is  clear  that  the  interaction  is  a  strong  one  and  the 
vortex  path  particularly,  when  the  vortex  is  near  the  blade  is  influenced  by  the  blade  and  the  interaction.  The 
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local  velocity  changes  significantly  during  the  interaction.  This  makes  techniques  based  upon  assumed  vortex 
position  or  shape  unlikely  to  provide  an  accurate  simulation. 

This  calculation  clearly  demonstrated  the  ability  of  the  code  in  accurately  convecting  a  vortex  that  was 
initially  located  far  upstream  of  the  blade,  downstream  toward  the  blade,  and  showed  the  flow  development  as 
the  vortex  interacted  with  the  blade. 


Figure  7:  Static  Pressure  contours  and  the  surface  pressure  coefficient  Cp  at  t  =  0.0 


Figure  8:  Static  Pressure  contours  and  the  surface  pressure  coefficient  Cp  at  t  =  4.0 


Figure  9:  Static  Pressure  contours  and  the  surface  pressure  coefficient  Cp  at  t  =  6.0 
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4.0  SUBMARINE  MAVUEVERS 

The  prediction  of  the  maneuvering  characteristics  of  underwater  vehicles  requires  the  accurate  calculation  of 
highly  complex  hydrodynamic  phenomena  which  are  not  amenable  to  simple  analytic  techniques.  Maneuvers 
of  most  bodies  are  difficult  to  predict  due  to  large  angles  of  attack  and  the  presence  of  unsteady  flow 
conditions.  Underwater  vehicles  pose  an  added  challenge  because  the  control  surfaces  used  for  generating  the 
necessary  forces  and  moments  during  the  maneuver  are  often  relatively  small  compared  to  the  overall  size  of 
the  buoyant  body  itself  Consequently,  viscous  forces  and  moments  generated  on  the  hull  and  appendages  can 
have  a  dominant  role  in  the  behavior  of  the  vehicle.  Additionally,  boundary  layer  and  vortex  wakes  generated 
by  the  hull  and  forward  appendages  interact  with  the  stem  appendages,  as  well  as  with  the  propulsor,  to 
provide  forces  and  moments  and  vehicle  dynamics  significantly  different  from  that  which  would  be 
experienced  without  such  interactions. 

The  propulsor  is  an  important  element  in  the  overall  maneuvering  character  of  the  vehicle,  and  it  adds 
considerable  further  complexity  to  any  predictive  effort.  Traditionally,  an  underwater  vehicle  propulsor  is 
located  at  the  stem  to  benefit  from  operating  in  the  vehicle's  wake.  However,  the  wake  region  introduces  a 
non-uniform  inflow  to  the  propulsor,  and  this  often  leads  to  out-of  plane  forces,  in  addition  to  the  direct  thmst 
and  moment  produced  by  the  propulsor.  These  out-of-plane  forces  can  produce  significant  moments,  because 
the  propulsor  is  often  located  a  considerable  distance  from  the  center  of  gravity  of  the  vehicle.  The  interaction  of 
flows  associated  with  the  hull,  appendages  and  propulsor  produces  a  highly  complicated  unsteady  hydrodynamic 
system. 

One  approach  for  prediction  of  the  trajectory  of  submarines  is  to  perform  numerical  simulations  based  on  the 
Unsteady  Reynolds- Averaged  Navier-Stokes  (UnRANS)  equations.  A  physics-based  prediction  of  the 
maneuvering  characteristics  can  be  obtained  by  coupling  forces  and  moments  data  from  these  solutions  with  a 
six  degree  of  freedom  (6-DOF)  rigid  body  motion  analysis.  Such  computer  simulations  can  provide  both 
quantitative  analysis  and  improved  understanding  of  flow  phenomena  affecting  both  design  and 
maneuverability.  This  information  can  contribute  to  both  improved  designs  and  the  safe  operation  of 
submarines.  Details  of  the  numerical  scheme  and  development  of  the  method  is  given  in  detail  in  [22].  Here 
validation  aspects  of  the  method  are  described. 

Flow  around  a  maneuvering  submarine  with  sailplanes,  stem  appendages,  and  a  rotating  propulsor  include 
many  vortical  stmctures  such  as  sail  vortex,  sailplane  vortices,  tip  vortices  emanating  form  stern  appendages 
and  propulsor  blades,  hull  vortex  feeding  sheet,  necklace  vortices  generated  in  the  root  of  the  sail,  and  in  the 
root  of  the  sailplanes.  These  vortices  affect  the  forces  and  moments  on  the  body  of  the  vehicle  and 
consequently  the  maneuvering  characteristics  of  the  vehicle.  Preservation  of  these  vortices  by  the  code  is 
therefore  essential.  A  vortex  preservation  computation  as  described  in  2.0  is  conducted  on  the  flow  solver 
code.  Based  on  that  calculation  appropriate  turbulence  model,  grid  density,  dissipation  factor,  and  time  step 
was  chosen  to  conduct  the  maneuvering  computations  [23]. 

4.1  VALIDATION  STUDIES 

Numerous  solutions  at  several  angles  of  drift  were  computed  for  a  submarine  configuration  (named  SUBOFF) 
with  four  stem  appendages.  Figure  10  shows  contours  of  the  x-component  of  vorticity  for  the  18°  drift  angle 
case.  Flow  traces  are  also  drawn  on  the  Figure  10.  The  x- vorticity  contours  together  with  the  flow  traces  show 
the  creation  of  two  strong  vortical  stmctures  emanating  from  the  hull.  The  computed  axial  force  and  lateral 
force  coefficients  are  also  compared  with  the  experimental  data  on  the  Figure  10.  The  computed  force  and 
moment  coefficients  are  in  extremely  good  agreement  with  the  measured  data  (much  better  than  anticipated  at 
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the  time),  even  up  to  the  maximum  drift  angle  of  18  degrees.  The  effects  of  the  capturing  vortical  structures 
are  clearly  indicated  in  the  axial  and  in  the  lateral  force  coefficients.  The  Baldwin-Lomax  turbulence  model  is 
used  in  these  calculations. 


Figure  10:  Contours  of  X-component  of  Vorticity  -  Experimental  Data  vs  Computational 


4.2  SIMULATION  OF  CRASHBACK  MANEUVER 

Deceleration  of  a  self-propelled  underwater  vehicle  by  means  of  reversing  the  angular  velocity  of  the  propeller 
is  called  a  crashback  maneuver.  This  maneuver  has  been  simulated  [22]  by  computing  a  starting  solution  for  a 
constant-speed  vehicle  and  then  decreasing  the  propeller  rotation  from  the  constant-rpm  forward  rotation 
speed,  through  zero  to  a  constant-rpm  reverse  rotation.  This  calculation  represents  the  most  severe  off-design 
condition  for  flow  through  the  propeller  and  demonstrates  capability  to  compute  a  very  complex  phenomenon 
regarding  both  the  flow  field  and  the  vehicle  dynamics.  During  this  maneuver,  the  vehicle  continues  to  move 
in  a  forward  direction.  Relative  to  the  vehicle  itself,  the  fluid  generally  flows  towards  the  stem,  with  path-lines 
passing  outboard  from  the  propeller.  The  reversed  rotation  of  the  blades  moves  the  fluid  near  the  propeller  in 
an  upstream  direction.  These  downstream  and  upstream  fluid  motions  create  a  shear  stress  and,  subsequently, 
generate  a  ring  vortex  located  just  outboard  and  downstream  of  the  propeller  blade  tips.  Figure  11  shows  this 
vortex  ring  at  some  instant  in  time.  Contours  of  the  flow  traces  vary  with  the  magnitude  of  the  velocity,  while 
the  blades  show  contours  representing  the  magnitude  of  the  static  pressure.  The  unsteady  motion  and  eventual 
decay  of  this  ring  vortex  generates  unsteady  side  forces  and  has  a  pronounced  effect  on  the  behavior  of  the 
vehicle.  The  time  history  of  propeller  speed,  vehicle  velocity,  forces,  and  circumferential  position  of  the  ring 
vortex  relative  to  the  blades  are  also  shown  in  Figure  11.  The  lateral  force  Fy  and  the  vertical  force  Fz  display 
large-amplitude,  large  wavelength  oscillations  between  positive  and  negative  values.  The  oscillations  are 
virtually  identical  for  the  two  force  components,  except  for  the  expected  90°  phase  difference.  These  out-of- 
plane  forces  also  have  small-amplitude  high-frequency  oscillations  (associated  with  the  blade-passing 
frequency)  superimposed  on  these  low-frequency  oscillations.  Figure  11  also  shows  a  trace  that  represents  the 
approximate  circumferential  position  of  the  low — pressure  region  of  the  ring  vortex  (or  rotating  stall  cell)  at 
selected  instants  in  time.  In  this  trace,  the  circular  symbols  indicate  the  approximate  circumferential  position 
of  this  low-pressure  region  or  cell,  measured  from  the  stem  plane  in  the  clockwise  direction  looking  upstream. 
Each  symbol  is  plotted  as  a  sine  of  its  angular  position  0,  and  the  dotted  line  is  a  least-squares  polynomial  fit. 
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Time 


Figure  11:  Vortex  Ring  Produced  by  the  Crashback  Maneuver  (ieft)  -  Time  History  of  Propeiier  Speed,  Vehicie 
Veiocity,  Forces,  and  Circumferentiai  Position  of  the  Ring  Vortex  Reiative  to  the  Biades  (right) 


4.3  SIMULATION  OF  DEPTH-CHANGING  MANEUVER 

Most  maneuvers  other  than  crashback  are  accompanied  by  movement  of  one  or  more  control  surfaces. 
Control  surface  movement  increases  grid  size  and  complexity,  due  to  control  surface  gaps  and  movement.  A 
simpler  but  approximate  approach  is  to  model  the  effect  of  control  surface  movement  by  imposing  external 
forces  to  fixed  control  surfaces.  This  approach  is  demonstrated  here  for  a  simple  depth-changing  maneuver  by 
externally  changing  the  magnitude  of  the  forces  and  moments  on  the  horizontal  stern  planes,  inducing  a 
pitching  moment  about  the  center  of  gravity,  and  this  leads  to  a  rise  or  dive  of  the  vehicle.  Figure  12  shows  the 
pressure  distribution  on  the  surface  of  the  fully  appended  SUBOFF  at  an  instant  of  time  during  the  rising 
maneuver  for  which  the  pitch  angle  is  about  23  degrees  and  the  angle  of  attack  is  about  13  degrees.  Since  the 
vehicle  is  pitching  about  its  center  of  gravity,  which  is  located  approximately  half  way  between  the  bow  and 
stern,  the  stagnation  point  on  the  nose  moves  toward  the  upper  part  of  the  nose,  and  the  stagnation  points  on 
the  tip  of  the  horizontal  stem  planes  move  toward  the  lower  surface  of  the  planes.  Once  the  vehicle  has 
elevated  to  a  certain  level,  the  applied  forces  on  the  horizontal  stem  planes  are  decreased  to  reverse  the 
direction  of  the  moment  from  a  pitch-up  moment  to  a  pitch-down  moment.  Figure  12  shows  the  pressure 
distribution  later  in  time,  after  the  vehicle  motion  has  reversed  from  pitch-up  to  a  pitch-down  course. 
Comparing  the  location  of  the  stagnation  points  on  the  nose  of  the  hull  and  the  tip  of  the  horizontal  stem 
planes  at  the  pitch-up  and  pitch  down  maneuvers,  it  can  be  seen  that  the  stagnation  points  have  shifted  to 
opposite  sides.  The  time  history  of  the  vehicle  inertial  position  and  the  applied  z-direction  force  on  the  top 
surface  of  the  stern  plane  are  shown  in  Figure  12. 
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Figure  12:  Surface  Pressure  Distributions  and  Trajectory  prediction  during  a  Depth-Changing  Maneuver 


5.0  SIMULATION  OF  A  REACTING  FLOW  IN  A  GASEOUS  TURBINE  COMBUSTOR 

Gas  turbine  combustors  are  designed  to  operate  with  stable  flames.  To  produce  a  good  mixing  of  the  reactants 
and  a  stable  flame,  the  reactants  usually  pass  through  some  swirlers  before  entering  the  combustion  chamber 
and  the  resultant  flow  in  the  chamber  is  a  very  high  swirling  flow.  The  simulation  presented  here 
demonstrates  the  capability  of  the  code  in  generating  and  maintaining  the  vortical  structures  inside  the 
combustion  chamber.  For  validation,  the  computed  results  are  compared  with  the  experimental  data. 

5.1  EXPERIMENTAL  SETUP  AND  THE  GEOMETRY 

The  schematic  of  the  experimental  model  gas  turbine  combustor  operating  on  air/methane  is  shown  in  Figure 
13.  The  experimental  work  is  performed  by  Bowman.  T.C.  and  Edwards,  C.  [24].  The  operating  condition  of 
the  gaseous  combustor  is  also  summarized  in  the  Figure  13.  The  overall  combustor  assembly  consists  of  three 
distinct  sections:  fuel  delivery,  main  combustion  chamber,  and  a  tailpipe.  Flow  is  delivered  through  two 
separate  co-annular  concentric  pipes.  The  low  velocity  methane  fuel  is  delivered  through  the  inner  pipe, 
whereas  the  higher  velocity  air  is  delivered  through  the  annulus  of  the  two  pipes.  Both  fuel  and  airflows  pass 
through  45°  helical  co-swirling  swirlers  and  become  highly  swirling  flows  as  they  enter  the  main  combustion 
chamber.  This  is  necessary  for  creating  a  lifted,  stable  flame  and  for  good  mixing  of  the  reactants.  A  close- 
up  of  the  45°  helical  co-swirling  swirlers  is  also  shown  in  Figure  13. 
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Figure  13:  Axial  cross  section  and  operating  condition  of  the  gaseous  combustor  operating  on  methane  fuel 


5.2  REACTING  FLOW  RESULTS 

The  numerical  simulations  performed  to  describe  the  flow,  solve  the  Reynolds-Averaged  Navier-Stokes 
equations  on  an  unstructured  discretization  with  added  turbulence,  chemical  species,  and  combustion  models. 
Turbulence  is  modelled  by  a  cubic  non-linear  k-epsilon  model  with  low  Reynolds  number  wall  integration 
[25].  The  chemistry-turbulence  interactions  are  represented  by  the  eddy-dissipation  combustion  model  of 
Magnussen  and  Hjertager  [26],  where  combustion  rate  is  assumed  to  be  controlled  by  the  turbulence  mixing 
rate  of  the  reactants. 

The  flow  structures  predicted  by  the  numerical  simulation  is  shown  if  Figure  14.  It  is  dominated  by  a  very 
strong  center  recirculation  zone,  and  by  another  ring  vortex  generated  in  the  upstream  comer  of  the 
combustor.  These  clearly  demonstrate  the  capability  of  the  code  in  predicting  these  flow  stmctures. 


Figure  14:  3-D  (left)  &  2-D  (right)  Flow  structures,  showing  Ring  Vortices,  Recirculation  zones 
Axial  Velocity,  and  particle  traces  (red  dashed  lines  are  vortex  cores) 

Figure  15,  shows  an  axial  cross  section  showing  the  computed  axial  velocity  contours,  azimuthal  contours, 
axial,  and  azimuthal  velocity  profiles  for  both  the  computed  and  the  measured  data  throughout  the  reacting 
field.  The  axial  velocity  is  high  in  the  air  stream  and  the  air  jet  quickly  spreads  toward  the  combustor  wall. 
The  central  recirculation  zone  around  the  centerline  with  a  length  approximately  equal  to  the  diameter  of  the 
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combustor  is  maintained.  Further,  another  recirculation  is  formed  close  to  the  wall.  These  recirculation  zones 
are  well  captured  by  the  simulation.  Re-circulating  hot  gases  in  the  center  and  the  bulkhead  corner  section 
causes  flame  stabilization. 


X[m] 

Figure  15:  Axial  &  Azimuthal  Velocity  Contours  &  Radial  Profiles 
Red  diamond  symbols  are  the  experimental  data,  solid  lines  are  the  computational  data, 
Pink  diamond  symbols  are  experimental  velocity  fluctuations  about  the  mean 
(They  are  used  as  the  50%  confidence  interval  about  the  mean) 


Good  agreements  between  the  computational  results  and  the  experimental  data  validate  the  suitability  of  the 
code  for  simulation  of  high  swirl  and  reacting  flows  in  combustion  chambers. 


6.0  SUMMARY  AND  CONCLUSION 

Any  computational  code  used  for  simulation  and  analysis  of  the  flow  in  or  around  real-world  devices  and 
vehicles  is  required  to  predict  the  flow  features  and  the  flow  quantities  associated  with  these  devices  and 
vehicles  with  reasonable  accuracy  and  in  a  reasonable  time  frame.  The  method  of  vortex  preservation 
presented  in  this  paper  offers  an  approach  where  the  time-accurate  computational  results  are  compared  with 
the  exact  solution  for  a  Lamb  type  vortex  convecting  in  a  free  stream.  Comparing  the  computed  results  with 
the  exact  solution  eliminates  the  uncertainties  that  exist  with  the  experimental  data  commonly  used  for  the 
validation  of  the  CFD  codes.  This  method  considers  effects  of  the  grid  density,  numerical  dissipation, 
turbulence  model,  and  grid  alignment  on  the  computed  results.  It  is  shown  that  if  a  numerical  scheme  can  not 
preserve  an  isolated  vortex,  it  can  not  be  used  for  investigation  of  the  flowfields  that  contain  these  flow 
structures. 

Three  real  world  scenarios,  each  with  a  different  flow  characteristics  are  presented  where  the  preservation  of 
the  vortical  structures  are  critical  to  the  outcome  of  the  simulations.  These  included  the  blade  vortex 
interaction  (BVI)  flow,  flow  around  self-propelled  maneuvering  submarines,  and  the  reacting  flow  in  a 
gaseous  combustor.  The  computed  results  were  compared  with  the  experimental  data  or  the  exact  solution, 
when  available.  In  many  cases  the  original  code  proved  to  be  incapable  of  preserving  the  flow  features 
experienced  in  the  real  world.  These  codes  were  accordingly  modified  to  include  these  capabilities.  The  final 
codes  used  for  the  simulations  were  developed  to  capture  and  maintain  the  flow  structures  associated  with 
each  case. 

A  new  method  was  validated  for  predicting  trajectory  of  appended  underwater  vehicles  with  rotating 
propulsors.  The  validation  work  included  trajectory  predictions  for  a  fully-appended  submarine  going  through 
a  crashback  maneuver  and  a  depth-changing  maneuver.  Furthermore,  the  results  of  a  number  of  computations 
performed  on  a  SUBOFF  geometry  with  four  stem  appendages  were  compared  with  experimental  data.  The 
comparisons  showed  excellent  agreements. 
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An  experimental  model  gas  turbine  combustor  exhibiting  the  typical  flow  features  associated  with  gas  turbine 
combustors  was  numerically  modelled  to  demonstrate  the  modelling  capabilities  of  another  computer  code  in 
predicting  the  flow  characteristics  of  these  combustors.  The  computed  results  were  validated  against  the 
experimental  data.  Flow  features  specific  to  gaseous  combustion  chambers  such  as  recirculation  zones  were 
also  captured  by  the  simulations. 
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